library(gdata)
setwd("/home/bruce/Dropbox/SenPE/Data/")
library(foreign)
sendat <- read.dta("senator_data.dta")

dat <- list()
files <- dir("./PEDat")

for(i in 1:length(files)){
	dat[[i]] <- read.xls(paste("./PEDat/",files[i],sep=""),stringsAsFactors=F, na.strings=c("NA",""))
}

# Names of senator columns in PE Data
scols <- paste("s",1:15,sep="")
# Congresses of dat
congs <- c(108,109,110,96+0:9)

# IDs for matching
id_master <- read.csv("id_master.csv",stringsAsFactors=F)


# Create the event Matrices

mats <- list()

nattend <- NULL
ind <- 1
ifow <- c(1,4:13)
for(i in 1:length(dat)){
	dati <- dat[[i]]
	congi <- congs[i]
	mati <- NULL
	if(is.element(i,ifow)){
	idsi <- id_master[which(id_master[,1]==congi),5]
	idsi <- idsi[which(!is.na(idsi))]
	}
	if(!is.element(i,ifow)){
		idsi <- as.character(unique(c(as.matrix(dati[,scols]))))
		idsi <- idsi[which(!is.na(idsi))]
		}
	veci <- numeric(length(idsi))
	for(j in 1:nrow(dati)){
		if(!is.na(dati$s1[j])){
			attend <- dati[j,scols]
			attend <- as.character(attend[which(!is.na(attend))])
			nattend <- c(nattend,length(attend))
			print(attend)
			vecij <- veci
			vecij[match(attend,idsi)] <- 1
			if(mean(vecij)!=0){
			mati <- rbind(mati,vecij)
			}
			}
		}
		colnames(mati) <- idsi
		mats[[ind]] <- mati
		ind <- ind+1
	}

wnets <- list()
# Create Weighted Networks
for(i in 4:13){
	mati <- mats[[i]]
	neti <- matrix(0,ncol(mati),ncol(mati))
	for(j in 1:nrow(neti)){
		for(k in (1:nrow(neti))[-j]){
			neti[j,k] <- sum(mati[,j]*mati[,k])
			}
		}
		colnames(neti) <- colnames(mati)
		rownames(neti) <- colnames(mati)
		wnets[[i-3]] <- neti
	}
	names(wnets) <- 96:105

library(igraph)
community <- list()

for(t in 1:length(mats)){
mati <- mats[[t]]
ntype <- c(rep(0,nrow(mati)),rep(1,ncol(mati)))
el <- NULL
for(i in 1:nrow(mati)){
	for(j in 1:ncol(mati)){
		if(mati[i,j]==1){
			n1 <- i-1
			n2 <- nrow(mati)+j-1
			el <- rbind(el,c(n1,n2))
			}
		}
	}

# Project the bipartite network

g <- graph.bipartite(ntype,c(t(el)))
g2 <- bipartite.projection(g)[[2]]
m <- as.matrix(g2)
edge_g2 <- cbind(m[4][[1]],m[3][[1]])
wg2 <- numeric(nrow(edge_g2))
for(i in 1:nrow(edge_g2)){
		indi <- edge_g2[i,1]
		indj <- edge_g2[i,2]
		wg2[i] <- sum(mati[,indi+1]*mati[,indj+1])
	} 

cg2 <- fastgreedy.community(g2)
cg2w <- fastgreedy.community(g2,weights=wg2)
mem <- community.to.membership(g2,cg2$merges,which.max(cg2$modularity))
memw <- community.to.membership(g2,cg2w$merges,which.max(cg2w$modularity))
community[[t]] <- list(cg2,cg2w,mem,memw)
}

fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")


comid <- numeric(nrow(id_master))

for(t in 4:length(mats)){
	congi <- congs[t]
	ind <- match(paste(colnames(mats[[t]]),congi,sep=""),paste(id_master$sen_pe,id_master$Congress,sep=""))
	comid[ind] <- community[[t]][[4]]$membership
	}
	
comid[which(id_master$Congress>105)] <- NA
	
write.csv(data.frame(id_master,comid),file="community_membership.csv")

## Get Modularity By Both Community and Party ##

mod_mat <- NULL

for(t in 4:length(mats)){
mati <- mats[[t]]
ntype <- c(rep(0,nrow(mati)),rep(1,ncol(mati)))
el <- NULL
for(i in 1:nrow(mati)){
	for(j in 1:ncol(mati)){
		if(mati[i,j]==1){
			n1 <- i-1
			n2 <- nrow(mati)+j-1
			el <- rbind(el,c(n1,n2))
			}
		}
	}
	

# Project the bipartite network

g <- graph.bipartite(ntype,c(t(el)))
g2 <- bipartite.projection(g)[[2]]
m <- as.matrix(g2)
edge_g2 <- cbind(m[4][[1]],m[3][[1]])
wg2 <- numeric(nrow(edge_g2))
for(i in 1:nrow(edge_g2)){
		indi <- edge_g2[i,1]
		indj <- edge_g2[i,2]
		wg2[i] <- sum(mati[,indi+1]*mati[,indj+1])
	} 

mem <- community[[t]][[4]]$membership
maxMod <- modularity(g2,mem,wg2)
memsi <- colnames(mats[[t]])
congi <- congs[t]
idmi <- subset(id_master,id_master$Congress==congi)
fowi <- subset(fowler,fowler$congress==congi)
fowid <- idmi$cosp_summary[match(memsi,idmi$sen_pe)]
party <- fowi$party[match(fowid,fowi$labels)]
party <- as.numeric(party==200)
party[which(is.na(party))] <- 0
partyMod <- modularity(g2,party,wg2)

mod_mat <- rbind(mod_mat,c(congi,maxMod,partyMod))

}

## Plot the degree distribution, with mean and median
detach(package:igraph)
for(i in 2:length(wnets)){
filei <- paste("/home/bruce/Dropbox/SenPE/Tex/degdist",names(wnets)[i],".pdf",sep="")
pdf(filei,height=4.5, width=5, family="Times", pointsize=13)
par(las=1,mar=c(6,6,.5,.5),cex.lab=2,cex.axis=1.75)
wni <- wnets[[i]] > 0
degi <- degree(wni,gmode="graph")
hist(degi,col="grey70",border="grey45",main="",xlab="",xlim=c(0,80),ylim=c(0,40),breaks=c(0:8)*10)
abline(v=c(mean(degi),median(degi)),lwd=2.75,lty=c(1,2))
title(xlab="Degree",line=3.5)
dev.off()
}
#Fowler's Data, Dem = 100, Rep = 200

fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")

# Plot communities based on party composition
library(colorRamps)
pdf("/home/bruce/Dropbox/SenPE/Tex/communities.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(96:105,seq(0,100,length=length(4:13)),col="white",xlab="",ylab="Cumulative Percent",bty="n",xlim=c(95.5,105.5),xaxt="n")
red2blues <- blue2red(100)[100:1]

for(i in 4:13){
	memsi <- colnames(mats[[i]])
	csizei <- community[[i]][[4]]$csize
	cmemi <- community[[i]][[4]]$membership
	coms <- 0:(max(which(csizei>1)-1))
	congi <- congs[i]
	csizei <- csizei[which(csizei>1)]
	props <- 100*csizei/sum(csizei)
	pdem <- numeric(length(coms))
	fowi <- subset(fowler,fowler$congress==congi)
	idmi <- subset(id_master,id_master$Congress==congi)
	for(j in 1:length(pdem)){
		cmemj <- memsi[which(cmemi==coms[j])]
		fowid <- idmi$cosp_summary[match(cmemj,idmi$sen_pe)]
		party <- fowi$party[match(fowid,fowi$labels)]
		pdem[j] <- mean(party==100,na.rm=T)
		}
	cszsrt <- csizei[order(pdem)]
	pdem <- pdem[order(pdem)]
	cols <- red2blues[1+round(pdem*99)]
	cssz <- cumsum(cszsrt/sum(cszsrt))
	lb <- c(0,cssz[1:(length(cssz)-1)])*100
	ub <- cssz*100
	for(j in 1:length(pdem)){
		rect(xleft=congi-.5,xright=congi+.5,ybottom=lb[j],ytop=ub[j],col=cols[j])
		}
	}

axis(1,at=96:105)
dev.off()

# Plot max modularity vs. polarization

pdf("/home/bruce/Dropbox/SenPE/Tex/modularity.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(mod_mat[,c(1,2)],ylim=c(min(c(mod_mat[,c(2,3)])),max(c(mod_mat[,c(2,3)]))), lwd=2.5,xlab="",ylab="Modularity",type="l")
lines(mod_mat[,c(1,3)],col="grey55",lwd=2.5)
abline(v=96:105,lwd=2.5,col="grey70",lty=3)
points(mod_mat[,c(1,2)],cex=1.5,pch=19)
points(mod_mat[,c(1,3)],pch=23,cex=1.5,col="grey55",bg="grey60")
legend("bottomright",legend=c("Max Modularity","Polarization"),col=c("black","grey55"),lwd=2.5,pch=c(19,23),pt.bg=c("black","grey55"),bg="white")

axis(1,at=96:105)
dev.off()

###### Comparison with Fowler ########

## Add centrality to Fowler's Data ##

fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")
fowid <- paste(fowler$labels,fowler$congress,sep="")

nc <- rep(0,nrow(fowler))
ec <- rep(0,nrow(fowler))
bc <- rep(0,nrow(fowler))
dc <- rep(0,nrow(fowler))

ncl <- rep(0,nrow(fowler))
ecl <- rep(0,nrow(fowler))
bcl <- rep(0,nrow(fowler))
dcl <- rep(0,nrow(fowler))

l_betweeness <- rep(NA,nrow(fowler))
l_connectedness <- rep(NA,nrow(fowler))
l_evcent <- rep(NA,nrow(fowler))
l_closeness <- rep(NA,nrow(fowler))
for(i in 1:length(wnets)){
library(igraph)
net <- wnets[[i]]
gnet <- graph.adjacency(net,mode="undirected",weighted=T)
congi <- as.numeric(names(wnets)[i])
distnet <- graph.adjacency(1/net,mode="undirected",weighted=T)
sp <- shortest.paths(distnet,algorithm="dijkstra")
csp <- c(sp)
csp_inf <- which(csp==Inf)
csp[csp_inf] <- 1.25*max(csp[-csp_inf])
msp <- matrix(csp,nrow(sp),nrow(sp)) 
new_cent <- 1/apply(msp,2,sum)
ev_cent <- evcent(gnet)
bet_cent <- betweenness(gnet)
cl_cent <- closeness(gnet)
pe_now <- colnames(net)
sm_now <- id_master[match(paste(pe_now,congi,sep=""),paste(id_master[,5],id_master[,1],sep="")),3]
sum_now <- id_master[match(paste(pe_now,congi,sep=""),paste(id_master[,5],id_master[,1],sep="")),2]
fow_now <- match(paste(sum_now,congi,sep=""),paste(fowler$labels,fowler$congress,sep=""))

nc[fow_now[which(!is.na(fow_now))]] <- new_cent[which(!is.na(fow_now))]
ec[fow_now[which(!is.na(fow_now))]]<- ev_cent$vector[which(!is.na(fow_now))]
bc[fow_now[which(!is.na(fow_now))]] <- bet_cent[which(!is.na(fow_now))]
dc[fow_now[which(!is.na(fow_now))]] <- cl_cent[which(!is.na(fow_now))]

fow_now[which(is.na(fow_now))] <- 1
betweeness <- fowler$betweeness[fow_now]
connectedness <- fowler$connectedness[fow_now]
evcent <- fowler$evcent[fow_now]
closeness <- fowler$closeness[fow_now]


if(i <length(wnets)){
	sum_next <- id_master[match(paste(sm_now,congi+1,sep=""),paste(id_master[,4],id_master[,1],sep="")),2]
	fow_next <- match(paste(sum_next,congi+1,sep=""),paste(fowler$labels,fowler$congress,sep=""))
	ncl[fow_next[which(!is.na(fow_next))]] <- new_cent[which(!is.na(fow_next))]
	ecl[fow_next[which(!is.na(fow_next))]]<- ev_cent$vector[which(!is.na(fow_next))]
	bcl[fow_next[which(!is.na(fow_next))]] <- bet_cent[which(!is.na(fow_next))]
	dcl[fow_next[which(!is.na(fow_next))]] <- cl_cent[which(!is.na(fow_next))]
	
	l_betweeness[fow_next[which(!is.na(fow_next))]] <- betweeness[which(!is.na(fow_next))]
	l_connectedness[fow_next[which(!is.na(fow_next))]] <- connectedness[which(!is.na(fow_next))]
	l_evcent[fow_next[which(!is.na(fow_next))]] <- evcent[which(!is.na(fow_next))]
	l_closeness[fow_next[which(!is.na(fow_next))]] <- closeness[which(!is.na(fow_next))]
	}

}

fowler <- data.frame(fowler,nc,ec,bc,dc,ncl,ecl,bcl,dcl,l_betweeness,l_connectedness,l_evcent,l_closeness)
fowler <- subset(fowler,is.element(fowler$congress,97:105))

# Scatter plot of centrality scores
for(i in 97:105){
filei <- paste("/Users/brucedesmarais/Desktop/Research/SenPE/Tex/twocents",i,".pdf",sep="")
nci <- subset(fowler$nc,fowler$congress==i)
mnci <- median(nci,na.rm=T)
connecti <- subset(fowler$connectedness,fowler$congress==i)
mconn <- median(connecti,na.rm=T)
cent_rho = cor(nci,connecti,use="complete.obs")
pdf(filei,height=4.5, width=5, family="Times", pointsize=13)
par(las=1,mar=c(6,6.25,.5,.6),cex.lab=2,cex.axis=1.75)
plot(connecti,nci,pch=4,col="grey35",xlab="",ylab="",cex=1.35)
abline(h=mnci,v=mconn,lwd=2)
title(xlab="Cosponsorship",line=3.5)
title(ylab="Press Events",line=4.75)
rr <- round(cent_rho,digits=3)
legend("topleft",legend=substitute(paste(rho," = ",rr),list(rr=rr)),col="white",cex=1.5,box.col="white",inset=-.05)
dev.off()
}

# Simple Means and Tests of Hypotheses
library(sciplot)
library(glm.nb)
med_nc <- median(fowler$nc)
med_ncl <- median(fowler$ncl)
ncid <- as.factor(fowler$nc > med_nc)
nclid <- as.factor(fowler$ncl > med_ncl)

filei <- "/Users/brucedesmarais/Desktop/Research/SenPE/Tex/barplot.pdf"
pdf(filei,height=4, width=5, family="Times", pointsize=13)
par(las=1,mar=c(4,4,2,.6),xpd=T)
bargraph.CI(ncid,fowler$pb,nclid,xaxt="n",legend=T,leg.lab=c(expression(paste("Low Centrality ", italic(t))),expression(paste("High Centrality ", italic(t)))),xpd=T,x.leg=1,y.leg=10.25,ylab="Mean Bills Passed",col=c("grey45","grey80"))
axis(1,at=c(2,5),labels=c(expression(paste("Low Cent. ", italic(t-1))),expression(paste("High Cent. ", italic(t-1)))))
dev.off()

setwd("/Users/brucedesmarais/Desktop/Research/SenPE/Data")
## Read in necessary libraries
# reads and writes Stata files
library(foreign)
# Implements negative binomial regression
library(MASS)

# Read in data
read.dta("fowler_plus.dta")

# Create an empty list in which to store model results
models <- list()

# Loop through Congresses, estimating a model for each
for(i in 97:105){
	# Estimate the model for the ith Congress
	modi <- glm.nb(pb~ncl*nc+connectedness+as.numeric(party==200)+seniority,data=subset(fowler,fowler$congress==i),x=T)
	# Store the ith model results
	models[[i-96]] <- modi
}
# Name the list by the corresponding Congress
names(models) <- 97:105

# write.dta(fowler,"fowler_plus.dta")

simple_models <- list()
for(i in 97:105){
modi <- glm.nb(cbind(pb)~nc,data=subset(fowler,fowler$congress==i),x=T)
simple_models[[i-96]] <- modi
}
names(models) <- 97:105

models_ni <- list()
for(i in 97:105){
modi <- glm.nb(cbind(pb)~ncl+connectedness+as.numeric(party==200)+seniority,data=subset(fowler,fowler$congress==i),x=T,maxit=50)
models_ni[[i-96]] <- modi
}
names(models) <- 97:105

models_nitm1 <- list()
for(i in 97:105){
modi <- glm.nb(cbind(pb)~nc+connectedness+as.numeric(party==200)+seniority,data=subset(fowler,fowler$congress==i),x=T,maxit=50)
models_nitm1[[i-96]] <- modi
}
names(models) <- 97:105


cvll.binom <- function(mod_obj,ydata){
	# Function to compute covariate contribution to CVLL
	x <- mod_obj$x
	y <- ydata[as.numeric(rownames(x)),]
	cvlls <- numeric(ncol(x))
	for(i in 1:(ncol(x))){
		xi <- x
		if(i>1){
			xi <- x[,-i]
			}
			cvli <- numeric(nrow(y))
			for(j in 1:nrow(y)){
				xt <- xi[-j,]
				yt <- y[-j,]
				xv <- xi[j,]
				yv <- y[j,]
				modij <- glm(yt~xt-1,family="binomial")
				lpij <- t(xv%*%cbind(coef(modij)))
				pij <- 1/(1+exp(-lpij))
				cvli[j] <- dbinom(yv[1],sum(yv),pij,log=T)
				}
				cvlls[i] <- sum(cvli)
		}
	cvlls
	
	}
	
cvll.pois <- function(mod_obj){
	# Function to compute covariate contribution to CVLL
	require(MASS)
	x <- mod_obj$x
	y <- mod_obj$y
	cvlls <- numeric(ncol(x))
	for(i in 1:(ncol(x))){
		xi <- x
		if(i>1){
			xi <- x[,-i]
			}
			cvli <- numeric(length(y))
			for(j in 1:length(y)){
				xt <- xi[-j,]
				yt <- y[-j]
				xv <- xi[j,]
				yv <- y[j]
				modij <- glm.nb(yt~xt-1,maxit=50)
				muij <- exp(xv%*%cbind(coef(modij)))
				cvli[j] <- sum(dnbinom(yv,size=modij$theta,mu=muij,log=T))
				}
				cvlls[i] <- sum(cvli)
		}
	cvlls
	
	}
	
cvl_contrib <- NULL
for(i in 1:length(models)){
	cvls <- cvll.pois(models[[i]])
	cvl_contrib <- rbind(cvl_contrib,cvls)
	}
	
	colnames(cvl_contrib) <- rownames(coef(models[[1]]))

cvl_contrib_ni <- NULL
for(i in 1:length(models)){
	cvls <- cvll.pois(models_ni[[i]])
	cvl_contrib_ni <- rbind(cvl_contrib_ni,cvls)
	}
	cvl_contrib <- cbind(cvl_contrib,cvl_contrib_ni[,1])
	
	
cvl_contrib_nitm1 <- NULL
for(i in 1:length(models)){
	cvls <- cvll.pois(models_nitm1[[i]])
	cvl_contrib_nitm1 <- rbind(cvl_contrib_nitm1,cvls)
	}	
	cvl_contrib <- cbind(cvl_contrib,cvl_contrib_nitm1[,1])
	


cvl_dif <- NULL
for(i in 2:9){
	difi <- cvl_contrib[,1]-cvl_contrib[,i]
	cvl_dif <- cbind(cvl_dif,difi)
	}

library(xtable)
rownames(cvl_dif) <- 97:105
colnames(cvl_dif)] <- colnames(t(coef(models[[5]])))[2:7]
xtable(t(cvl_dif))

# Ropeladder plots of coefficients
# Horizontal barplots of cvll contribution
# CI plots of interaction effects


# Effect Simulation 

eff_sim <- function(mod){
	require(mvtnorm)
	x <- mod$x
	clast <- x[,2]
	cnow <- x[,3]
	xmu <- apply(x[,-c(1,2,3,7)],2,mean)
	xmu <- matrix(rep(xmu,100),nrow=100,ncol=length(xmu),byrow=T)
	qlast <- quantile(clast,prob=c(.1,.9))
	snow <- seq(min(clast),max(clast),length=100)
	betaMu <- coef(mod)[c(2,3,7)]
	betaCov <- summary(mod)$cov.unscaled[c(2,3,7),c(2,3,7)]
	betas <- rmvnorm(1000,mean=betaMu,sigma=betaCov)
	#betas <- cbind(1,betas[,c(1,2)],coef(mod)[4],coef(mod)[5],coef(mod)[6],betas[,3])
	lplq <- NULL
	lpuq <- NULL
	xlq <- cbind(1,qlast[1],snow,xmu,qlast[1]*snow)
	xuq <- cbind(1,qlast[2],snow,xmu,qlast[2]*snow)
	for(i in 1:1000){
		betai <- betas[i,]
		lpli <- (betai[2]+betai[3]*snow)*sd(cnow)
		lplq <- cbind(lplq,lpli)
		}
		return(list(snow,100*(exp(lplq)-1),clast))
	}

eff_sims <- list()
for(i in 1:length(models)){
	eff_sims[[i]] <- eff_sim(models[[i]])
	}

eff_dists <- list()
for(i in 1:length(eff_sims)){
	qlq <- t(apply(eff_sims[[i]][[2]],1,quantile,prob=c(.025,.975)))
	qlq <- cbind(apply(eff_sims[[i]][[2]],1,mean),qlq)
	#quq <- t(apply(eff_sims[[i]][[3]],1,quantile,prob=c(.025,.975)))
	#quq <- cbind(apply(eff_sims[[i]][[3]],1,mean),quq)
	eff_dists[[i]] <- qlq
	}

names(eff_dists) <- 97:105
for(i in 1:length(eff_dists)){
filei <- paste("/Users/brucedesmarais/Desktop/Research/SenPE/Tex/effects",names(eff_dists)[i],".pdf",sep="")
pdf(filei,height=4.5, width=5, family="Times", pointsize=13)
par(las=1,mar=c(6,6,.5,.5),cex.lab=2,cex.axis=1.75)
lq <- eff_dists[[i]][,2]
uq <- eff_dists[[i]][,3]
plot(eff_sims[[i]][[1]],eff_dists[[i]][,1],col="black",lwd=4.5,type="l",ylab="",xlab="",ylim=c(min(lq),max(uq)))
abline(h=0,lwd=5,col="grey55")
lines(eff_sims[[i]][[1]],eff_dists[[i]][,1],col="black",lwd=4.5,lty=1)
lines(eff_sims[[i]][[1]],eff_dists[[i]][,2],col="black",lwd=4.5,lty=2)
lines(eff_sims[[i]][[1]],eff_dists[[i]][,3],col="black",lwd=4.5,lty=2)
#lines(eff_sims[[i]][[1]],eff_dists[[i]][[2]][,1],lwd=4.5)
#lines(eff_sims[[i]][[1]],eff_dists[[i]][[2]][,2],lwd=4.5,lty=2)
#lines(eff_sims[[i]][[1]],eff_dists[[i]][[2]][,3],lwd=4.5,lty=2)
rug(eff_sims[[i]][[3]],ticksize=.07,col="grey50")
title(xlab="Lagged Centrality in PE",line=3.5)
title(ylab="Percentage Change",line=4.5)
dev.off()
}

# Cosponsorship Community Detection
# Ropeladders of effects
# Table of Predictive Performances


##### Cosponsorship Community Detection #####

## Make the lists of cosponsorship matrices, sponsorship counts, and member ids

senmats <- list()
senspons <- list()
senmems <- list()
for(i in 93:108){
	if(i < 100){
		matfile <- paste("/home/bruce/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_matrices/senate_matrices/0",i,"_senmatrix.txt",sep="")
		idfile <- paste("/home/bruce/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_members/senate_members/0",i,"_senators.txt",sep="") 
	}
	if(i >= 100){
		matfile <- paste("/home/bruce/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_matrices/senate_matrices/",i,"_senmatrix.txt",sep="")
		idfile <- paste("/home/bruce/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_members/senate_members/",i,"_senators.txt",sep="")
	}
	senmati <- read.csv(matfile,header=F)
	senmati <- (senmati==1) + (senmati == 2)
	senmems[[i-92]] <- read.csv(idfile,header=F)
	amat <- matrix(0,nrow(senmati),nrow(senmati))
	sponn <- numeric(nrow(senmati))
	for(j in 1:nrow(amat)){
		for(k in (1:nrow(amat))[-j]){
			amat[k,j] <- sum(senmati[j,]*senmati[k,])
			amat[j,k] <- sum(senmati[j,]*senmati[k,])
			
			}
		}
		print(i)
	senmats[[i-92]] <- amat
	senspons[[i-92]] <- sponn
}


# Make list of node-level data

sendat <- read.csv("/home/bruce/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/SH.csv",header=T,stringsAsFactors=F)
sendats <- list()
for(i in 93:108){
	sendati <- subset(sendat,is.element(sendat$congress,i))
	sendats[[i-92]] <- sendati[match(senmems[[i-92]][,3],sendati$ids),]
}

### Community Detection ###

cosp_comm <- list()

for(t in 4:13){
g2 <- graph.adjacency(senmats[[t]],weighted=T,mode="undirected")
cg2w <- spinglass.community(g2)
memw <- cg2w$membership
cosp_comm[[t]] <- list(cg2w,memw)
}

com_mod <- NULL

for(t in 4:13){
g <- graph.adjacency(senmats[[t]],weighted=T,mode="undirected")
mem <- cosp_comm[[t]][[1]]$membership
maxMod <- modularity(g,mem,weight=E(g)$weight)
senmemi <- as.character(senmems[[t]][,1])
congi <- congs[t]
idmi <- subset(id_master,id_master$Congress==congi)
sum_idi <- idmi$cosp_summary[match(sub("  "," ",senmemi),idmi$cosp_senmem)]
fowi <- subset(fowler,fowler$congress==congi)
party <- sendats[[t]]$party
party <- as.numeric(party==200)
party[which(is.na(party))] <- 0
partyMod <- modularity(g,party,weight=E(g)$weight)
com_mod <- rbind(com_mod,c(congi,maxMod,partyMod))

}

# plot modularity vs. polarization in cosponsorship

pdf("/Users/brucedesmarais/Desktop/Research/SenPE/Tex/com_mod.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(com_mod[,c(1,2)],ylim=c(min(c(com_mod[,c(2,3)])),max(c(com_mod[,c(2,3)]))), lwd=2.5,xlab="",ylab="Modularity",type="l",yaxt="n")
lines(com_mod[,c(1,3)],col="grey55",lwd=2.5)
abline(v=96:105,lwd=2.5,col="grey70",lty=3)
points(com_mod[,c(1,2)],cex=1.5,pch=19)
points(com_mod[,c(1,3)],pch=23,cex=1.5,col="grey55",bg="grey60")
legend("bottomright",legend=c("Max Modularity","Polarization"),col=c("black","grey55"),lwd=2.5,pch=c(19,23),pt.bg=c("black","grey55"),bg="white")
axis(1,at=96:105)
axis(2,at=(3:10)/100)
dev.off()

# Plot colored community membership for cosponsorship

library(colorRamps)
pdf("/Users/brucedesmarais/Desktop/Research/SenPE/Tex/cosp_comm.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(96:105,seq(0,100,length=length(4:13)),col="white",xlab="",ylab="Cumulative Percent",bty="n",xlim=c(95.5,105.5),xaxt="n")
red2blues <- blue2red(100)[100:1]

for(i in 4:13){
	memsi <- colnames(mats[[i]])
	csizei <- cosp_comm[[i]][[1]]$csize
	cmemi <- cosp_comm[[i]][[1]]$membership
	coms <- which(csizei>1)-1
	congi <- congs[i]
	csizei <- csizei[which(csizei>1)]
	props <- 100*csizei/sum(csizei)
	pdem <- numeric(length(coms))
	for(j in 1:length(pdem)){
		cmemj <- memsi[which(cmemi==coms[j])]
		party <- sendats[[i]]$party[which(cmemi==coms[j])]
		pdem[j] <- mean(party==100,na.rm=T)
		}
	cszsrt <- csizei[order(pdem)]
	pdem <- pdem[order(pdem)]
	cols <- red2blues[1+round(pdem*99)]
	cssz <- cumsum(cszsrt/sum(cszsrt))
	lb <- c(0,cssz[1:(length(cssz)-1)])*100
	ub <- cssz*100
	for(j in 1:length(pdem)){
		rect(xleft=congi-.5,xright=congi+.5,ybottom=lb[j],ytop=ub[j],col=cols[j])
		}
	}

axis(1,at=96:105)
dev.off()

### Ropeladder Plots of Effects ###
effnames <- c("lpecent","pecent","cospcent","party","senior","peinter")
for(i in 2:7){
	t <- i-1
	filei <- paste("/Users/brucedesmarais/Desktop/Research/SenPE/Tex/",effnames[t],".pdf",sep="") 
	pdf(filei,height=3.15, width=2.75, family="Times", pointsize=13.5)
par(las=1,mar=c(2.5,2.75,.5,1),cex.lab=1.25)
	effmat <- NULL
	for(j in 1:length(models)){
		effj <- summary(models[[j]])$coef[i,1]
		sej <- summary(models[[j]])$coef[i,2]
		lowj <- effj-1.96*sej
		hij <- effj+1.96*sej
		effmat <- rbind(effmat,c(effj,lowj,hij))
		}
		cols <- paste("grey",50*as.numeric(sign(effmat[,2])!=sign(effmat[,3])),sep="")
		plot(effmat[,1],1:9,ylim=c(.75,9.25),xlim=c(min(effmat[,2]),max(effmat[,3])),pch=19,col=cols,cex=1.25,ylab="",xlab="",yaxt="n")
		abline(v=0,lwd=1.5,lty=2)
		abline(h=1:9,col="grey50",lty=3,lwd=2.5)
		points(effmat[,1],1:9,pch=19,cex=1.25,col=cols)
		segments(x0=effmat[,2],x1=effmat[,3],y0=1:9,y1=1:9,lwd=3.5,col=cols)
		axis(2,at=1:9,labels=97:105)
		#title(xlab="Coefficient",line=2.25)
		dev.off()
	}

